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ABSTRACT 






An error ellipse plotting routine was developed to 
provide real-time indication of Kalman filter performance. 

The study included an evaluation of the Hewlett-Packard HP-86 
computer system's capability for providing real-time 
tracking information and an evaluation of the computer's 
possible use on the three-dimensional underwater tracking 
range at the Naval Underwater Weapons Engineering Station, 
Keyport, Washington. A series of tracking runs were used 
to demonstrate both linear and extended Kalman filtering. 
Information obtained from the error ellipses was used to 
modify filter parameters for improved filter perf orm'ance . 

It was found that the error ellipse was useful as a tool 
for indicating filter performance and for making decisions 
regarding filter parameter modification. The HP-86 provided 
accurate, reliable results and it could be used for on-line 
graphics. However, the computing speed fo the HP-86 computer 
as used in this study was too slow for on-line processing 
of the three-dimensional tracking problem. 
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I. INTRODUCTION 



The Kalman filter's importance as an estimaror and 
predictor is well documented. Providing real-time informa- 
tion concerning filter performance so that on-line adjust- 
ments to filter parameters can be made continues to be an 
area of high interest. This study investigates the 
usefulness of the error ellipse as a tool for providing a 
real-time indication of filter performance. 

Part of the investigation involves an evaluation of the 
Hewlett-Packard HP-86 computer system's capacity to operate 
in a real time tracking environment, and its capabilities for 
providing information concerning filter performance. The 
rationale behind this investigation is based on the require- 
ment for the Naval Underwater Weapons Engineering Station, 
Keyport, Washington, to accurately acoustically track 
torpedoes on a three-dimensional underwater tracking range. 

A knowledge of the range operation is not within the scope 
of this study. It is sufficient to know that presently 
the range receives four time measurements every 1.31 seconds, 
and these measurements are nonlinear functions of the 
torpedo position. 

To gain a better understanding of the error ellipse, a 
4-state tracking scenario was chosen for this study. 
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Initially, the linear tracking problem is discussed, 
followed by an investigation of the nonlinear problem. Of 
primary interest are on-line methods to improve filter 
performance using information provided by the error ellipse 
for filter parameter modification. 
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II. KALMN FILTER THEORY 



A. LINEAR MATHEMATICAL MODEL 

1. The Plant 

For this model the state and measurement equations 
for the plant are linear. Hence the discrete form is used. 
The assumed plant model is described by a linear, vector 
difference equation: 

x(k+l) = _^(k) + ^(k) + rw(k) (State Equation) (2-1) 
and a linear, vector measurement equation: 

^(k) = cx(k) +v(k) (2-2) 

where : 

x(k) is an n-dimensional column vector, denoting the 
state of the plant at "time" k. 

u(k) is the deterministic control input, an m-vector, 
at time k. 

w(k) is a p-dimensional vector representing any random 
forcing inputs at time k. 

^(k) is a q-dimensional vector representing measure- 
ments made at time k. 

v(k) is a q-dimensional vector representing random 
measurement made at time k. 

r_, and c are assumed constant coefficient matrices of 
dimension nxn, nxm, nxp , and qxn respectively. 
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Noise Processes 



In order to place probabilistic structure on the 
noise processes v(k) and w(k) the following assumptions are 
made : 

(a) v(k) and w(k) are individually white processes, 
that is, for any k and 1, with k 1, v(k) and 
v(l) are independent random variables, and w(k) 
and w(l) are independent random variables. 

(b) v(k) and w(k) are individually zero mean, 
Gaussian processes with known covariances. 

(c) v(k) and w(k) are independent processes. 

Thus for the measurement noise: 



MEAN : 



E[v(k)] = 0 (k = 0,l,2,3,. . .) 



(2-3) 



COVARIANCE: E[v(k) v'^( 1) ] = E[v(k) ]E[v^(l) ] 



A 



?k 



k 7^ 1 

k = 1 



or E[v(k)v^ (1)] = (k,l=0,l,2, . . . ) 



(2-4) 



where 




is the Kronecker delta function defined as: 
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kl 



0, k 7^ 1 

1, k = 1 
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Likewise for the random forcing input: 



ECw(k)] =0 (k=0 ,1,2 ,3 , . . . ) (2-5) 

ECwCk)^"^ (1)] = 5^^ k,l = 0 ,1,2 ,3 , . . . ) (2-6) 

Q and R are nonnegative definite symmetric for all 
k. Also since v(k) and w(k) are zero mean and independent 
then : 



E[v(k)w'^ (1)] = 0 



(2-7) 



For the purposes of this study, unless otherwise 

specified Q and R, are considered to be known and constant, 

although both may be time varying . 

3 . Initial State Description 

For the initial state of the difference equation (2-1) 

it is unlikely that ^ will be available. Hence, it is 

assumed that x is a Gaussian random variable of known mean 

— o 

X and known covariance P , i.e., 

— o — o 



E[x(0)] = X 

— — o 



E{ [ (x 
— o 



X ) ] [ (x 
— o — o 



X )]^} 

— o 



= p 
— o 
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This choice for the initial state has the advanrage 
of causing the subsequent estimation scheme to be unbiased 
for all t. [Ref. l] Further it is assumed that the initial 
state and the measurement noise are uncorrelated: 

E[x(0)v"^ (k)] = E [v( k)x"^ ( 0 ) ] = £ ( k= 0 , 1 , 2 , 3 , . . . ) 

Also the initial state and the random forcing input are 
uncorrelated : 

E[x(0)w'^ (k)] = E[w(k)x'^ (0)] = 0 (k=0 , 1 , 2 , 3 , . . . ) 

B. DISCRETE-TIME ESTIMATION 

1 . The Estimator Equations 

The estimation problem involves generating an 
optimal estimate for x ( j ) for the system described by the 
difference equation (2-1) from the noisy measurements 
0 ),£( 1 j ) . This estimate will be denoted by 

A 

x(j/j), which means the estimate of x at time j given 
measurements at times up to and including time j . The 
estimate must be optimal in the sense that the expected 
value of the sum of the squares of the error in the estimate 
is a minimum , i , e . : 

E{[(x(k/k) - x(k)]'^[x(k/k) - x(k)]} = minimum 
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The estimator is characterized by the linear relationship: 



x(k/k) = x(k/k-l) + G(k)[z(k) - cx(k/k-l)] (k=0,l,2,..) 



( 2 - 8 ) 



where 

A 

x(k/k) is the optimal (minimum variance) estimate of 

x(k) given observations at times up to and 
including k. 

A 

x(k/k-l) is the optimal one-step prediction of x(k) 
given observations at times up to and 
including k-1. 



G(k) is the optimal estimation gain matrix which 

will minimize the variance of estimation 
error . 



For the initial estimate x(0/0), the estimator 

A 

equation (2-8) is initialized with x(0/-l), which is not a 
random variable. If x(0/-l) is selected such that: 



x(0/-l) = E[x(0)] = X 



it can be shown that this choice of x(0/-l) makes the 
estimator unbiased for all k. [Ref. 1] The estimators best 
available information concerning x(k-l) is the estimate 

A 

x(k-l/k-l) , therefore it is reasonable to assume that 

x(k/k-l) = _^(k-l/k-l) + ^(k-1) (2-9) 

is the best prediction. 



13 



In summary, equations (2-8) and (2-9) are the estimator 

A 

equations, v;ith x(0/~l) " initial condition. 

2 . Gain and Covariance Equations 

Without going into detailed derivations, the optimal 
estimator gains, G(k) , used in the estimator equation (2-8), 
are those which satisfy: 

G(k) = P(k/k-l)c'^[^(k/k-l)C^ + R]"^ (2-10) 



P(k/k) = [I - G(k)C]P(k/k-l) 
P(k+l/k) = _^(k/k)({)^ + Q 



with the initial conditions: 



P(0/-1) = P 

Q 



E{ [x(0) 



X ] [x ( 0 ) 

— o — 






where 



P(k/k) 



E{[x(k/k) - x(k)][x(k/k) - x(k)]"^} 



( 2 - 11 ) 

( 2 - 12 ) 



is the covariance of estimation error matrix. 
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P(k/k-l) = E{[x(k/k-l) - x(k)][x(k/k-l) - x(k)]‘} 

is the covariance of one-step prediction error matrix. 

^ = E[^(k) .w(k) (k)-r^ (k)] 

is the state excitation matrix. 

P(k/k) and P(k/k-l) are symmetric, positive definite matrices. 

Several observations can be made concerning the linear 
Kalman gain (2-10), covariance (2-11, 2-12) and estimator 
(2-8) equations. 

(a) The estimator gains, G(k), do not depend on the 
measurement data and hence can be precomputed ,* stored , and 
used as the processing measurements become available. 

(b) Although not obvious from the equation, the 
time-varying gain, G(k) , depends in time as: 

G(k) = [Ref. 2] (2-13) 

Thus the effect is to weight the correction term, 

[£(k) - cx(k/k-l)], in the estimator equation (2-8) less 
heavily as time progresses. The advantage of a greater 
initial weight allows for possibly large differences 
between £(k) and x(k/k-l) during the initial observations, 
and a large gain will result in a significant change in the 
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next estimate. This advantage is also borne out in that 
there is less confidence in the quality of the estimates 
during the early observations compared with the quality 
after numerous observations. Hence the later an observation, 
the less drastic an estimate will be altered or affected by 
an isolated observation discrepancy. 

(c) In general, the variance of estimation error 
decreases in a manner analogous to the gain schedule (2-13), 
i.e., it decreases as k grown larger, reflecting greater 
confidence in the estimate as the number of observations 
increases. Selection of the proper initial condition, , 
is important when studying the effect of measurement errors 
on the behavior of the estimate. So P should be assigned 
pessimistic values which would correspond to a lack of 
information about the initial state. In cases when the 
initial state is completely unknown, then P^->“»I. [Ref. 3] 

(d) The Q matrix serves to compensate for model 
errors and prevents the covariance matrix from becoming too 
small or optimistic. A small covariance matrix would result 
in a small filter gain, and subsequent observations are 
essentially ignored, which could result in filter divergence. 
The Q_ matrix prevents G(k) from approaching zero by adding 
uncertainty to the system which is reflected in a 
degradation of certainty (increase in P(k+l/k)). 
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C. NONLINEAR ESTIMATION - EXTENDED KALxMAN FILTER 

In many practical applications, the state equations 
and/or measurement equations are nonlinear. Before the 
Kalman filter equations can be used, the problem must be 
linearized and the Kalman filter equations are applied with 
some modification. 

1 . Nonlinear Model 

Consider a nonlinear discrete system of state and 
observation equations given by: 

x(k+l) = f(x(k), u(k) ,k) + w(k) (2-14) 



and 



^(k) = h(x(k) , k) + v(k) (2-15) 

In these equations f and h are nonlinear functions 
of the state variable x, w(k) is the plant excitation noise, 
and v(k) is the measurement noise. The plant noise and 
measurement noise are assumed to be uncorrelated, zero-mean, 
and white. The same equations (2-3 thru 2-7) apply as for 
the linear model. 

2 . Extended Kalman Filter Equations 

In order to apply the linear filter equations, 
equations (2-14) and (2-15) are expanded about the best 
estimate of the state at that time and only the first-order 
terms are kept. 
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That is 5 defining A(k) as: 
(x(k/k) , u (k ) 5 k) 



(x(k/k-l ) ) 

seen from the above equations , the filter 
and x(k/k-l) are used as the "best" 
estimates about which the linearization is performed. The 
matrices A(k) and H(k) must be used to generate G(k) so it 
is available to process £(k) when it is obtained. The 
modified extended Kalman filter equations are then: 

Gain Equation: 

G(k) = P(k/k-l)H^ (k)[H(k) 

Filter Update Equation: 

x(k/k) = x(k/k-l) + G(k)[z(k) - h(x(k/k-l))] (2-17) 



• P(k/k-l) • (k) + R]"^ 

(2-16) 



A(k) = ^ 



8f 

3x 



and 



H(k) 



9h 

9x 



As can be 

As 

estimates, x(k/k) 
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Prediction Equation: 



x(k+l/k) = f(x(k/k), u(k) , k) 



Covariance of Estimation Error Equations: 



P(k/k-l) = A(k-l)P(k-l/k-l)A (k-1) + Q(k-l) 
P(k/k) = [I - G(k)H(k) ]P(k/k-l) 



(2-18) 

(2-19) 

( 2 - 20 ) 



/V 

For the initial estimate x^O/O)? equation (2-17) is 
initialized with x(0/-l) with 



x(0/-l) = E[x(0)] 



X 

— o 



x(0/-l) is also used to initially evaluate H(k). 



As in the linear case: 
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III. ERROR ELLIPSOIDS 



A. THEORY 

Since the estimate x(k/k) is unbiased in the Kalman 

filter equations, the P(k/k) matrix represents the covariance 

of the error in the estimate. If the estimate were biased, 

P(k/k) would represent the second-moment matrix rather than 

the covariance matrix. Hence, P(k/k) provides significant 

information about the accuracy of the estimate. If the 

physical model is accurately described by the state and 

measurement equations (2-1, 2-2), then P(k/k) can be used 

to describe the manner in which the estimate converges 

(or diverges) to the true state. Examination of the P(k/k) 

matrix directly, element by element, is not a realistic 

2 

approach, since the matrix contains n elements, where n is 
the number of state variables . To simplify the situation 
the concept of the error ellipsoid is used. [Ref. 4] 

As discussed earlier, the assumptions are made that the 
initial state of the plant is Gaussian, as are the random 
processes v(k) and w(k) . Using these assumptions it follows 

A 

that x(k) and x(k/k) are also Gaussian since they are linear 
combinations of Gaussian variables and deterministic 
quantities. Using the same rationale, the estimation error, 
defined as: 
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/V 



e(k/k) 



= x(k/k) 



X (k ) 



is also Gaussian. Using the fact that the mean of the 
estimation error is 0, the probability density function for 
£(k/k) is: 

p^[e(k/k)] = P(k/k) I 

exp[-l/2e'*' (k/k)P ^(k/k) e (k/k) ] (3-1) 

The density function, p^[e(k/k)] will have a constant 
value whenever the exponent has a constant value. That is: 

-l/2e'^ (k/k)P ^ (k/k)e(k/k) = c 



or 



e^ (k/k)P ^(k/k)e(k/k) = c^, (3-2) 

where c is an arbitrary constant. 

As demonstrated by Sorenson [Ref. 1] and Kirk [Ref. 2], 
it can be shown that the locus of points e(k/k) which 
satisfy equation (3-2) are hyperellipsoids. For the 
two-dimensional case which is of concern, equation (3-2) 
describes an ellipse. This can be seen by fixing time and 
rewriting (3-2) as: 
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T 2 

e we = c 



(3-3) 



where 



.-1 



w = P (k/k) (a 2 X 2 symmetric matrix) 



Expanding the left side of (3-3) gives 



^ 22®2 ‘ ^ 



which because of symmetry gives 



^ 11®1 ^ 22®2 ^ 



(3-4) 



Since w^^> 0 , W22>0? ^11 ^22^^12’ equation (3-4) 

describes an ellipse, in which the principal axes do not 
coincide with the coordinate axes . The ellipse is rewritten 
in terms of y(l) and y(2) as the coordinate axis, and it can 
be shown that y(l) and y ( 2 ) are the eigenvectors of w with 
and defined as the corresponding eigenvalues. [Ref. 2] 
(See Figure 3-1). The equation for the ellipse can be 
rewritten in terms of a coordinate system having unit vectors 
in the directions of y(l) and y ( 2 ) as the basis vectors. 

The ellipse equation becomes: 

X^y^(l) + X^y^(2) = c^ (3-5) 
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Figure 3-1 Error Ellipse 



Remembering that w = P ^(k/k) , it can be shown that the 
corresponding eigenvectors and eigenvalues for w ^=P(k/k) 

'2 X 



are y(l), y(2), a,, and , where a. = ^ and a„ = 



Equation (3-5) can be rewritten: 



y ^ ( 1) y ^ ( 2 ) ^ 

+ = c (3-6) 

“l “2 
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For the purposes of this study, the term error ellipsoid 
refers to the specific case when c = 1. So equation (3-5) 
becomes : 

y ^ ( 1) y ^ ( 2 ) 

+ = 1 

a 2 

In terms of the Cartesian coordinates x’, y', which use e^ 
and e 2 a-s basis vectors: 

x’ = XCOS0 + ysin9 
y’ = -xsinS + ycos9 

where 9 is the angle of rotation of the axes and can be 
computed from: 

I 2cov ( e^ e ) 

9 = 1/2 tan~’'^[ ] [Ref, 5] 

var(e ) - var(e ) 

X y 



For a given value of c it is possible to integrate the 
probability density over the surface of the error ellipse to 
obtain the probability that a particular sample point will 
lie within the ellipsoid. For this study, n = 2, c = 1, 
and the probability the error is inside the ellipse is 0.394. 

In summary, the error ellipsoid can be used to 
characterize the concentration of the estimate about the 
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true value of the state. A decrease in the magnitude of an 
axis of the ellipse is an indication that the error in the 
estimate is decreasing in that direction. 

One important item that needs to be pointed out is that 
often the components of the state vector represent entirely 
different types of variables, for example the components 
might represent range, velocity, and depth. Since the two- 
dimensional ellipses are determined by using two components 
of the state vector, it is reasonable to examine submatrices 
relating state variables of the same character. Doing so 
will preclude most scaling difficulties when plotting the 
ellipses, and provide more meaningful insight in to the 
results . 

B. ERROR ELLIPSOIDS AND FILTER DIVERGENCE 

Thus far the discussion has centered around using 
submatrices of the P(k/k) covariance of estimation error 
matrix as the input for the error ellipse to indicate 
filter performance. As proposed by Heffes [Ref, 6] and 
Nishimura [Ref. 7], P(k/k) can be considered as a "design" 
covariance matrix p'^ , using the assumption that the only 
errors are in Q^, R^, and P^ with the following inequalities 
holding for all k: 



, R^ > R^ 

— k — ^k ’ — k k 



-o — — o 



(3-7) 
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with the subscript ’’d” indicating designed and "a" indicating 
actual. Equation (3-7) implies more input noise, more 
measurement noise, and more initial state uncertainty in the 
design than actually exists. This conservative filter design 
results in a somewhat pessimistic design error covariance 
P*^(k/k) . The actual error covariance P'^(k/k) resulting 
from using a filter designed with Q , and P^ is related 
to the design error covariance in the following manner: 



P^(k/k) > P^(k/k) 



This result is particularly useful when one simply does 
not know accurately the noise covariance of the input or 
output, but an upper bound is known. Designing assuming 
the noise covariance is at its upper bound will result in 
P^(k/k) being upper bounded by P'^Ck/k) . In some sense a 
worst case design results. Filter divergence exists when 
the design error covariance P'^Ck/k) remains bounded while 
the error performance matrix P^(k/k) becomes very large 
relative to p'^Ck/k) or is, in fact, unbounded. 
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PROBLEM DEFINITION 



A. PROBLEM DESIGN 



The purpose of the tracking problem is to study the use 
of error ellipsoids as real-time indicators of filter 
performance. In order to keep the design model realistic 
albeit reasonably simplified for ease of study, a 
two-dimensional tracking problem using several different 
tracks has been selected. 

1 . Linear Tracking 

All tracks are based on an x-y coordinate system 
with the target moving in the x or y direction relative to 
the sensor located at the origin. Thus for aircraft 
tracking, altitude is considered constant, as is depth for 
torpedo tracking. 

Defining the state variables as: 
x^ = X x-coordinate of the target location. 

X 2 = X velocity of target in x-direction. 

x^ = y y-coordinate of the target location. 



x^^ = y velocity of target y-direction. 
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resulting in a state .vector : 



X 



X 



« 

X 

y 

y_ 



( 4 - 1 ) 



The following are the state equations: 
x^(t) = X2(t) 

x^Ct) = w^Ct) 

x^Ct) = x^(t) 

x^(t) = W2(t) 



( 4 - 2 ) 



where w^(t) and W2(t) are assumed to be uncorrelated, 
random processes that account for unknown target accelerations 
and nonlinear target motions. Writing the discrete form of 
the state equations gives: 
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I 



T ^ 

x^(k+l) = x^(k) + T • x^Ck) + Y~ ’ 



X 2 (k+ 1 ) = x^(k) 




T 


• w^(k) 










x„(k+l) = x^(k) 

O V-/ 




T 


* x^(k) + — • w^(k) 


X(^(k+1) = x^(k) 


+ 


T 


• w^(k) 



or 



x(k+l) = ^(k) + rw(k) 



(4-3) 



(4-4) 



With T, the sampling period equal to 1 second, the matrices 
c}) and r are : 





“l 


1 


0 


o“ 




“ .5 


0“ 




0 


1 


0 


0 


r — 


1 


0 




0 


0 


1 


1 


i - 


0 


.5 




0 


0 


0 


1 




0 


1 



It is assumed that the sensor gives noisy, but 
uncorrelated measurements of x and y. Hence the discrete 
measurement equations are: 
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z^(k) = x^(k) + v^(k) 
z„(k) = x-(k) + v„(k) 



( 4-6) 



with Vj^(k) and v^Ck) uncorrelated random noise. 
Thus for the measurement equation: 



z(k) = cx(k) + v(k) 



(4-7) 



the matrix c is: 



c 



10 0 0 



0 0 10 



For the initial run and unless otherwise noted the 
following values for the Gaussian random processes will be 
used : 



E[v(k) ] = 0 for all k 



E[v(k)v (k)] = 



20x10 0 



20x10 



m = R for all k 



30 



a 

— V 



150 



150 



m = the standard deviation of 
measurement noise. 



E[w(k)] = 0 for all k 



E[w(k)w (k)] = 



100 



10 0 



2 2 

(m/sec ) = cov w for all k 



a 

— w 



10 



10 



m/sec = the standard deviations of the 
random forcing input. 



The covariance of estimation error matrix is initialized: 



P =P(0/-1) 
— o — 



0 

0 

0 



0 0 0 

10 ^ 0 0 

0 10 ^ 0 

0 0 10 ^^ 
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Since the filter is to be unbiased, the initialization: 
x(0/-l) = 21 q " Initial condition of the problem. 

2 . Nonlinear Tracking 

The state equations are the same as for the linear 
tracking problem. The measurement equation is considered 
as a noisy range measurement by the tracking sensor and is 
characterized as: 

z(k) = [x^(k) + x^(k)]^^^ + v-|(k) (4-8) 



Thus z(k) is a nonlinear function of the states. Using 
equation (2-14) and with T = 1 second: 



f (x(k) , u(k) , k) 



x^(k) + X 2 (h) 
x^(k) 

X 2 ( k ) + X ( k ) 
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Taking partial derivatives of f with respect to x gives: 



8f 



A(k) = 



3x 



(x(k/k ) 5 u( k ) 5 k) 



110 0 
0 10 0 
0 0 11 
0 0 0 0 






using equation (2-15): 



h(x(k) , k) = Cx^(k) + Xp(k)]^'^^ 



and taking the partial derivative with respect to x gives 



9h 



(x(k/k-D) 



x^ ( k) 



0, 



27/7727777172 

x^ (k) +x^ (k ) J 



x^Ck) 

[x^(k)+X3(k) 




x(k/k-l) 
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Using the above results with the same values for the 
random noise processes as previously stated for the linear 
case, the extended Kalman filter equations can now be applied 
to the nonlinear tracking problem. 

B. COMPUTER SIMULATION 

1 . Computer Hardware/Software 
a . Hardware 

All computer simulations were run on the Hewlett- 
Packard HP-86 personal computer. This particular model was 
chosen to evaluate its capabilities in determining its 
usefulness in actual torpedo tracking at the underwater 
tracking range at Naval Underwater Weapons Engineering 
Station, Keyport , Washington. The HP-86 system used included 
keyboard, 9 inch CRT monitor connected through an integrated 
monitor interface, and two HP Flexible Disc Drives connected 
through an integrated disc interface. Plotting was done on 
a HP-7225B Plotter and printing on a HP-2631G Printer. These 
peripherals were interfaced using a HP-IB Interface module. 
Because the system uses interface select codes, the HP-13 
factory preset code was set at 7 , which is the select code 
for the printer/disc interface. This code however did not 
work when interfacing with the external plotter and printer, 
since duplicate select codes are not allowed. So the 
internally set select code of the HP-IB was set to 8 for 
proper system operation. 
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b. Software 



The HP-86 has 60K built-in, useable bytes of 
computer memory, expandable to 572K using either 32K, 64K, 
or 128K Memory Modules. A HP-86 plug-in ROM was required 
to operate the external plotter. Also a Matrix ROM was used 
to reduce program length and computer run time . 

All programs were written in BASIC programming 
language using REAL (full) precision, which provides 15 digit 
precision. Appendix B provides an explanation of the program 
options and Appendix C contains the program listings used for 
this study. 

2 . Track Generation 

To evaluate the real-time use of the error ellipse as 
an indicator of filter performance, Monte Carlo simulation 
runs were made. Four tracks were generated by separate 
programs and one second incremental values of x, y, v^ , and 
v^ were stored in data files . Appendix I contains an 
explanation of the generation of tracks three and four. 

3 . Noise Generation 

In order to simulate the random noise processes, the 
computer’s random number generator was used and the generated 
numbers scaled accordingly. For each track and each 
different value of noise sigma, a different generator "seed 
number" was used. These noise values produced were added to 
the applicable true track values to simulate a sensor 
measurement corrupted by independent Gaussian noise. For all 
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cases where filter parameters were varied for a particular 
track under a specific noise condition, one noisy track was 
generated, stored, and used throughout that particular 
simulation. This was done for ease of filter performance 
comparison . 

4 . Gating Scheme 

In order to preclude catastrophic filter failure due 
to excessive measurement noise, a bound was established for 
the maximum acceptable limits of measurement noise. A 
three-sigma gate was designed using the covariance of the 
measurement noise, R, and the predicted covariance of error 
matrix P(k/k-l) . The gate is defined as: 

Gate(k) = 3(p.. (k/k-1) + 

^11 11 
max 

This gate is the maximum error allowable for the measurement 
at time k. If the absolute difference between the actual 
measurement received and the predicted measurement is 
greater than the three-sigma gate, then that particular 
measurement data is rejected as unacceptable. When this 
occurs, the filter gain, G(k) , is set to zero, resulting in 
that measurement being ignored and the prediction of the 
states set equal to the estimate, that is: 

x(k/k) = x(k/k+l) 
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5. 



Collection of Statistics 



In order to study error ellipses as an indicaror of 
filter performance, statistics were calculated, on line, 
after each measurement during the Monte Carlo run. The 
statistics computed were relative error (in some cases the 
error was normalized), error mean, error variance, and error 
covariance for the positional variables. The following 
equations apply: 



Relative Error = e(k/k) = x(k) - x(k/k) 



(filter error residual) 



k 

Error Mean = e(k/k) - 1/k Z e^(j/j) 

j=l 



Error Variance = VarCe^(k/k)] 



1/k Z [e . ( j / j ) ]-[e . (k/k) ] 
. , 1 1 



i=x 



1 ’ 




X 



4 



Positional 

Error 

Covariance 

Matrix 



VarCe (k/k)] 

X-, 



k 

1/k Z [e (j/j)][e (j/j)] 

:=i ^1 ^3 

-[.e (k/k)*e (k/k)] 

X, x^ 



k 

1/k Z [e (j/j)][e (j/j)] 

j=l ^1 ^3 

-[e (k/k)*e (k/k)] 

X-, x„ 



VarCe (k/k)] 

Xo 
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V. TARGET TRACKING AND ERROR ELLIPSE ANALYSIS 



A. LINEAR TRACKING 

Track 1 depicts a target approaching at a constant 

velocity of 223.6 feet per second. The solid line of Figure 

5.1 indicates the true track of the target and the numbers 

along the track indicate the time in seconds. The target 

was tracked in a measurement noise, a = 150, with the 

’ -7 ’ 

random forcing noise 0 =1. R was set for 20,000. The 

— w — 

dots indicate the filtered track using the linear Kalman 

filter equations. Figures 5.2 and 5.3 are the filter error 

ellipses for this run, computed at increments of 10 seconds. 

As can be seen the ellipse size decreases with increasing 

time indicating filter convergence. The computed ellipse 

surface areas shown on the figures confirm this. Figure 5.4 

shows the filtered track for the case where a has been 

—V 

increased to 300 and all other parameters remain the same. 

The error ellipses of Figure 5.5 computed for 10-second 
increments show increasing area indicating filter divergence. 
Figure 5.6 shows the error ellipses for the same track run 
but this time the ellipses were computed using a "statistics 
window" of 10. By this is meant that the ellipses were 
derived from statistics computed for the last 10 data 
measurements. All previous data is disregarded. Using this 
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method, the ellipses of Figure 5.6 show filter convergence 
from iterations 15 to 25. The filtered track of Figure 
5.4 confirms rhis . Figure 5.7 shows the resulrs for rhe 
same track parameters, except in this case a statistics 
window of 5 was used. The window 5 ellipse area at time 
25 (5179 sq ft) is much less than the area of either the 
run with the statistic window of 10 (11,540 sq ft) or the 
run with no window at all (65,500 sq ft). This is expected 
since the filter is essentially "locked on" at time 20, 
and the window 5 ellipse at time 25 disregards all data 
previous to time 20. Figure 5.8 shows the error ellipses 
for the same track but the measurement noise was increased 
to ^ = 400, while R was kept at 20,000. The error ellipses 
indicate filter divergence, and indeed the filtered track 
headed off in the wrong direction. 

Track 2 depicts a target approaching at a constant speed 
of 500 feet per second in the -y direction. Figure 5.9 
depicts the solid line track and the dots indicate the linear 
filtered track with cr = 150, a =1 and R = 20,000 . Figure 

— V — 

5.10 and 5.11 are the error ellipses for the run. No 
statistics window was used. Other than the fact that the 
ellipse area is decreasing, the shape of the ellipse provides 
little additional information. Figure 5.12 and 5,13 show the 
ellipses for the normalized error. With the same measurement 
noise sigma for both the x-position and y-position measurements, 
the normalized error ellipse's shape and orientation reflect 
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the target track proximity to either axis. In Figure 5.12 
the ellipse major axis indicates a large normalized error in 
the x-direction. This is to be expected since the target 
maintains a constant x-position of 500 feet, while the 
y-position is initially 20,000 feet. However, near time 40 
as the target approaches the x-axis , the normalized y error 
increases and becomes so large at the x-axis crossing that 
the normalized x error becomes insignificant in comparison. 
This can be seen in Figure 5.13 for the ellipses at time 
45 and 55 compared with the ellipse at time 35. The 
normalized error ellipse is an excellent indicator of 
target proximity to an axis , but the rapid shifts in ellipse 
surface area make it difficult to determine filter 
convergence . 

Track 3 depicts a target approaching on a parabolic 
track at a speed of 200 feet per second. Figure 5.14 is the 
true track, with a linearly filtered track indicated by the 
dots. For this run o =150, a =10, and R=20,000. Figure 5.15 
are the error ellipse plots for times 40, 50, and 60, the 
period of the highest rate of change in x- and y-velocity. 

The ellipse areas increase with time indicating divergence. 
Figure 5,16 and 5.17 are the ellipse, plots using a 10-data 
point and a 5-data point statistics window respectively. 

In both cases the ellipse areas for time 60 are less than 
time 50 indicating the filter has tracked around the curve. 
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Using error ellipses based on a statistics window in this 
tracking situation provides a better indicator of filter 
performance. 

For the second run of track 3, was increased to 300 
and all other parameters remained the same. Figure 5.18 
shows the resulting filtered track, which obviously didn'i: 
track around the curve. With the large amount of measurement 
noise and considering that the maximum random forcing input, 
i.e. the maximum acceleration in the x- and y-direction, 
occurs between times 40 and 60, it is a logical place to lose 
track. Another factor to be considered is the decrease in 
gain as k increases. By time 40 the gains have little 
influence. So a system was incorporated in the program to 
"reinitialize” the filter by setting the gains to G(0), if 
certain conditions were met. After several trial and error 
runs, it was determined that if a statistics window of 10 
were used, and the gains reinitialized if the error ellipse 
area increased consecutively a certain number of times, 
that the filtered track would follow around the curve. 

Figure 5.19 is the filtered track using a- statistics window 
of 10 and reinitializing the filter if the error ellipse area 
increased consecutively during five 1-second increments. 

The error ellipses for that run are shown in Figures 5.20 
and 5.21. As indicated on the figures, the filter was 
reinitialized four times and the ellipse areas for the 
10-second increments increased, until reinitializing the 



41 



filter at time 52 locked the filter in. Consequently the 
error ellipse areas for time 60 and 70 decreased, indicating 
convergence. 

In the final run, the filter was reinitialized after 10 
consecutive 1-second error ellipse increases, with all other 
parameters remaining the same. Figure 5.22 is the filtered 
track; Figures 5.23 and 5.24 are the error ellipse plots 
for this run. As indicated on Figure 5.24 the filter was 
reinitialized only once at time 57 . A comparison of ellipse 
areas for the last two runs shows that, with the exception 
of time 70, the areas were larger in the first run when the 
filter was reinitialized 4 times. This is expected since 
reinitializing results in larger gains producing more widely 
vary estimates, and hence greater error variance. At time 
70 the error ellipse area of the first run is less since in 
the first run the last filter reinitialization occurred at 
time 52 versus time 57 in the second run. The first 
filtered track had more time to settle out by time 70, 
resulting in less error. 

B. NONLINEAR TRACKING 

Track 4 depicts a target moving at a constant velocity 
of 50 feet per second (30 kts) in the x-direction for 15 
seconds, at which time the target turns and travels in the 
-y-direction . (See Figure 5.25) Using the extended Kalman 
filter with £_^^ = 30 and R = 900 , a series of tracking runs were 



2 

made for various values of COVW (a =from 20 to 200). In none 

w 

of these runs did the filter successfully track the target 
around the turn. Figure 5.26 shows the results for the case 
when COVW=150. In this instance the filter lost track as the 
target came out of the turn. Trying to track through a turn 
using a constant COVW (and hence, a constant Q) did not work. 
So a scheme was devised to vary COVW dependent on information 
derived from the error ellipse. After several trial runs for 
this particular track, it was determined that if the error 
ellipse area increased consecutively for 7 iterations of 
k, COVW would be doubled, and if the ellipse area decreased 
consecutively for 5 iterations of k, COVW would be halved. 
Initially the the trial runs were made without a statistics 
window, and the filter did not track successfully. Without 
the statistics window, the old data weighted down the statis- 
tics, and the error ellipses were not indicative of what was 
currently happening. So it was decided to use a statistics 
window. Windows of 5, 10, and 15 were tried. Window 5 
proved to be too responsive and window 15 not responsive 
enough. So a statistics window of 10 was chosen for the 
tracking run. With P^ = 10^, £_^-30, R = 900 , and COVW initially 
set at 20, the tracking run was made. Figure 5.27 depicts 
the filtered track output, and Figures 5.28-5.30 are the 
10-second incremental error ellipses for the run. As 
can be seen, the filter did track around the curve. Figure 
5.29 shows the ellipse areas are becoming less between times 
55 and 65. Also indicated below the plots are the values of 
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COVW for the k time the plot was computed. During this 
particular run COVW varied from an initial value of 20 up 
to 160, and then decreased to 40 by the end of the run. 

Using the same filter parameters as above except was 
increased to 200, and R to 40,000, another run was made. The 
filter did not track at all. During this run COVW varied 
from an initial value of 20 up to 40 and decreased to 5 by 
the end of the run. Obviously, the criteria for increasing 
and decreasing COVW was not effective. More trial runs were 
necessary to determine the optimum consecutive increases 
or decreases of the ellipse areas before adjusting COVW 
accordingly . 

A second approach to the nonlinear tracking problem, one 

that was used earlier for the linear case, is to reinitialize 

the filter if certain conditions are met. Again, using 

trail and error runs with and without statistics windows, 

it was determined that using a statistics window of 10 gave 

the best results. Using initial conditions of COVW-150, 

2 

P =10 , a =30, and R=900, several runs were made, reinitializ- 

— o —V — 

ing the filter if the error ellipse area increased consecu- 
tively for a certain number of iterations of k. Of the runs 
attempted, the best results were obtained when the filter was 
reinitialized if the area increased for 5 consecutive 
iterations of k. Figure 5.31 is the filtered track for this 
run, and Figures 5 . 32-5 . 34 are the error ellipse plots. 
Indicated below the plots are the times, k, when the filter 
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was reinitialized. For this particular run the filter was 

reinitialized 5 times. It should be noted that v;hen 

2 

reinitializing Tihe filter 5 ?(k/k-l) was reset to IG x 

was not large enough to be effective in getting the 

2 

filter back on track. The initial value of 10 was used for 

P to reflect the high confidence in the initial state 
— o ^ 

conditions. Other values of P did not work as well. 

— o 

Using the same parameters as above except the filter was 
reinitialized after 7 consecutive error ellipse area 
increases, another run was made with the resulting track 
depicted in Figure 5.35. A comparison with Figure 5.31 
reveals that using 7 area increases as the criterion for 
filter reinitialization resulted in poorer filter performance, 
as can be shown from the error ellipses. 

The final nonlinear filter tracking runs attempted 
involved simultaneously varying COVW and filter reinitializa- 
tion. The results were disastrous, and highly unpredictable. 
It was an interesting experiment in futility, and no 
meaningful results could be obtained. 
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Figure 5.1 Solid Line Track 1, Vel 223.6 ft/sec; Dots 

Indicate Filtered Track, a =150, a =1, R=20,000 

— V — — 



46 



155 



Y-POS (FT) 




COVW=l SIGV=150 R=20000 TIME AREACSQ FT) LEG 

10 1057E+001 

20 8790E+000 

30 7540E+000 



Figure 5.2 Filtered Track 1 Error Ellipses at 10 Second 
Increments, a =150 
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Figure 5.3 Filtered Track 1 Error Ellipses at 10 Second 
Increments, a =150 
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Figure 5.4 Filtered Track 1, a =300, a =1, R=20,000 
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Figure 5.5 Filtered Track 1 Error Ellipses for a_ 
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Figure 5.6 Filtered Track 1 Error Ellipses for o^=300, 
Using Statistics Window of 10 
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Filtered Track 1 Error Ellipses for ^ =300, 
Using Statistics Window of 5 ^ 



5 2 



Figure 5,7 
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Figure 5.8 Filtered Track 1 Error Ellipse for a =400 
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Figure 5.9 Solid Line Track 2, Vel 500 ft/sec, Dots Indicate 

Filtered Track for =150, o_ =1, R=20,000 
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Figure 5.10 Filtered Track 2 Error Ellipses, a^=150 
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Figure 5.11 Filtered Track 2 Error Ellipses, o_^^=150 
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Figure 5.12 



Filtered Track 2 Error Ellipses, =150, Error 
Normalized 
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Figure 5.13 



Filtere(i Track 2 Error Ellipses, o_^ = 150, Error 
Normalized 
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Figure 5,14 Solid Line Track 3, Vel 200 ft/sec, Dots 

Indicate Filtered Track for o_ =150, a =10, 
R= 2 0 ,00 0 ^ 
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Figure 5.15 Filtered Track 3 Error Ellipses, ^ =150 
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Figure 5.16 Filtered Track 3 Error Ellipses, a^^ = 150, 

Statistics Window 10 
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Figure 5.17 Filtered Track 3 Error Ellipses, o_ =150, 

Statistics Window 5 ^ 
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Figure 5.18 Filtered Track 3 for a =300, a =10, R=20,000 
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Figure 5.19 Filtered Track 3 for a =300, a =10, R=20,000, 
° — V — w — 

Statistics Window 10, Reinitialized at 5 

Consecutive Ellipse Area Increments 
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Figure 5,20 Filtered Track 3 Error Ellipses, o_^=300 , 

Statistics Window 10 , Reinitialize 5 
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Figure 5.21 Filtered Track 3 Error Ellipses, a_ =300, 

Statistics Window 10, Reinitialize 5 
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Figure 5.22 Filtered Track 3 for o =300, a =10, R=20,000 
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Statistics Window 10, Reinitialize 10 
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igure 5.23 Filtered Track 3 Error Ellipses, a^=300, 
Statistics Window 10, Reinitialize 10 
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Figure 5.24 Filtered Track 3 Error Ellipses, a. =300, 
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Statistics Window 10, Reinitialize 10 
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Figure 5 .,25 Track 4, Vel 50 ft/sec 
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Figure 5.26 Filtered Track 4 
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Figure 5.27 Filtered Track 4, a =30 , R=900, Statistics 
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Figure 5.28 Filtered Track 4 Error Ellipses, Statistics 

Window 10, Varying COVW 
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Figure 5.29 Filtered Track 4 Error Ellipses, Statistics 

Window 10, Varying COVW 
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Figure 5.30 Filtered Track 4 Error Ellipses, Statistics 

Window 10, Varying COVW 
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Figure 5.31 



Filtered Track 4, a =30, R=900, a^=150, 
Reinitialize 5 



76 



6000 




SIGV=30 R=900 COVW=150 


TIME 


AREA (SO FT) LEG. 


RESET 




15 


1805E-001 


0 




25 


1211E+000 


23 




35 


5743E+000 


33 



Figure 5.32 Filtered Track 4 Error Ellipses, a_ =30, 

Reinitialize 5 
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Figure 5.33 Filtered Track 4 Error Ellipses, £.^ = 30, 

Reinitialize 5 
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Figure 5.34 Filtered Track 4 Error Ellipses 
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VI. CONCLUSIONS AND ^RECOMMENDATIONS 



A. ERROR ELLIPSE 

The filter error ellipse proved useful as a tool for 
indicating filter performance. The information provided by 
the ellipse, particularly surface area changes, was used to 
make decisions concerning the alteration of the filter 
parameters. Several approaches for using the error ellipse 
were applied to both the linear and nonlinear tracking 
problem and the results are summarized below: 



Procedure Applicable Comments 

Fil ter 

( Linear or 
' Extended) 



Statistics Both 
Window 



Useful in keeping the error 
ellipse current and responsive 
to present data. The ellipse 
reflects most recent data; old 
data is disregarded. Beneficial 
when making a decision concerning 
filter parameter modification. 
Normally a better indicator of 
filter convergence or divergence 
than without a statistics window. 



Normalized Both 

Error 

Ellipse 



Aid in displaying error trends 
as target approaches coordinate 
axis or origin. Not practical 
in determining filter conver- 
gence/divergence due to rapid 
ellipse changes in vicinity of 
axes . 
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Reinitialize Linear- 

Filter- Set G(0) 

Increasing Extended-Set 

Ellipse Area PCk+l/k)= 

10^‘ X P 
— o 



Use increasing ellipse size as an 
indicator of divergence, and as a 
decision-making device to 
reinitialize filter. Parricularly 
valuable later in track when gains 
and P(k+1/K) have settled out. 

More effective when used in 
conjunction with statistics window. 



Error Extended 

Ellipse 

Expansion/ 

Compression 
to vary COVW 
(Adaptive Q) 



This technique of varying COVW 
based on error ellipse area in- 
crease or decrease is particularly 
useful in a tracking environment 
containing large variations in the 
random forcing input. The proce- 
dure is more effective when used 
in conjunction with a statistics 
window. Using this technique 
with filter reinitialization was 
unsuccessful . 



B. COMPUTER PERFORMANCE 

The HP-86 proved to be an extremely reliable computer 
with no downtime experienced during the 5 month period of 
operation. Full (Real) precision was used throughout the 
study providing 15 digit precision, which was more than 
adequate. The use of the Matrix ROM reduced program length 
by 90% and increased computing speed by a factor of about 6. 
Although not used in this study, a Statistics ROM would have 
undoubtedly further increased computing speed. For any 
further study using the HP-86, it is recommended that a 
Statistics ROM be procured. 

It took approximately 2 seconds for each incremental 
time measurement data to be sequenced through the filter 
equations, both for the linear and nonlinear case. This 
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computing time also included all statistics computations, 
with an incoming data rate of 1 set of measurement data per 
second, this computing speed is not sufficient for on-line 
processing. As previously mentioned, the Naval Underwater 
Tracking Range receives a series of 4 measurement times 
sequentially every 1.31 seconds. The range's three- 
dimensional tracking problem will necessarily involve more 
than the 4 state variables used in this study. Hence 
greater matrix dimensions resulting in longer computing 
times can be expected. 

The HP-86 CRT graphics were used extensively to provide 
error ellipse plots during the tracking runs. Using a "no 
frills" approach to plotting, i.e. plotting without x-y 
axis or labelling, it took approximately 2.5 seconds per 
ellipse plot. The ellipse plotting routine used involved 
sines and cosines, plotted point by point in 30 degree 
increments for a total of 360 degrees. This method was 
somewhat slow. Had there been available a graphics program 
that would sketch in the ellipse around the intersected 
major and minor axis, the graphics presentation could have 
been speeded up. But since ellipse plotting for every 3 to 
5 increments of time provided sufficient "real-time" 
information, the time of 2.5 seconds per plot was tolerable. 

Summarizing, the HP-86 could be used to compute 
statistics and provide graphics in the real-time underwater 
tracking environment, if the graphics were required not more 
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< 



often than 3 to 5 seconds. However, before the HP-86 can be 
considered feasible for real-time Kalman filter processing, 
more investigation is needed in finding ways to speed up 
computer processing time such as parallel processing, 
additional use of manufacturer-provided ROMs and machine 
language programming. 
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APPENDIX A 



TRACK GENERATION 



1. TRACK THREE 

Target movement is in the x-y plane v/ith the tracking 
sensor located at the origin of the cartesian axes. The 
target follows a parabolic track (see Figure A-1) at a 
constant speed of 200 feet per second. 




Figure A-1 Track 3 
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The parabolic equation is: 



= 4p(x-h) 



where p = 1000 and h = 1000, resulting in: 
= 4000(x-1000) 



Initial target location is (x,y) = (8000, 5291.5) 
x-direction velocity is given by: 

v^ = V cos(Angle) 

and target y-direction velocity is: 



v^ = V sin(Angle) 



where the Angle is obtained 



Angle = tan"^(|^) 



and V is obtained from: 



, 2 ^ 2 , 1/2 

V = (v + V ) 

X y 



from : 



. Target 



(A-1) 



(A-2) 
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where the argument of the inverse tangent is the slope of a 
small increment (less than 1 second) at each successive data 
point. Using a sampling period of one second, the data 
points (x,y) can be obtained from: 



x(k+l) = x(k) + 


v.^ ( k ) 


y(k+l) = y(k) + 


V (k) 

y 



Table A-1 gives the numerical values for the four states. 



2 . TRACK FOUR 




Target movement 


simulates a 30-knot torpedo (velocity 50 


feet per second) at 


a constant depth. The target's initial 


position is (x,y) = 


(3250, 5000), and it is moving in the v 


direction with v =0, 

y 


. (See Figure A-2.) The target remains 



on a straight course for 15 seconds, and then executes a 90 
degree turn and travels in the -y direction at v = -50 ft/sec 

. y 

The trajectory of the target turn is described by a 90 degree 
arc of a circle of radius, R = 1000, with the circle centered 
at (x,y) = (4000, 4000). The 90 deg arc will be traversed 
in : 

2ttR/4v sec. = IOtt sec. (where v = 50 ft/sec) 
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21 
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30 

31 

32 

33 
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35 

36 

37 

38 

39 

40 
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42 

43 

44 
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TABLE A-1 



NUMERICAL VALUES OF STATES FOR TRACK THREE 
X X-VEL Y Y-VEL 



8000. CO 
7813. 10 
7626.17 

7439.58 
7253.33 
7067.44 
6881.93 
6696. 62 

6512.14 
6327.90 

6144.14 
5960. 87 

5778. 14 
5595.98 
5414.4 1 

5233.49 
5C53.25 

4873. 74 
4695.02 
4517. 15 
4340. 19 
4164. 22 
3589.31 
3815. 57 
3643.09 
3472. 00 
3302.43 

3134.54 
2968. 5C 
2804. 52 
2642.85 
2483. 76 

2327. 59 

2174. 74 
2C25.66 
1880. 95 
1741.28 

1607. 50 
1480.67 
1362. 08 
1253.38 

1156. 65 
1C74.64 
1011. 20 
1000. CO 



-186.90 

- 186. 93 
-186.60 

- 186. 25 
-185. 89 
-185.5 1 
-185 . 1 1 

- 184. 68 

- 1 84 . 2 4 

- 183. 76 
-183 . 26 

- 182. 73 

- 182. 17 

- 181. 57 
-180.92 

- 180. 24 
-179.51 

- 178. 72 

- 177 . 87 
-176. 96 
-175. 97 

- 174. 9 1 

- 173 .74 

- 172. 48 
-171 .09 

- 169. 57 
-167 .89 
-166. 04 
-163. 98 

- 161. 67 
-159. 09 
-156. 17 
-152. 86 

- 149 . 07 
-144 .72 

- 139. 67 

- 133 .78 

- 126. 8 3 
-118.59 
-108.70 

-96 . 73 
-82. 0 1 
-63 . 44 
-37. 24 
0.00 



5291 . 50 
5220. 38 
5148. 27 
5075* 26 
5001. 33 
4 9 26 . 4 3 

4850. 54 

4773. 60 

4695. 59 
4616. 45 
4536. 14 

4454. 60 
4371 . 79 

4287. 65 
4202. 10 
4115. 09 

4026. 54 
3936. 37 

3844. 49 
3750. 81 
3655. 24 

3557. 65 
3457. 93 
3355. S3 

3251. 52 

3144. 52 
3034. 75 
2922. 01 
2806. 06 

2686. 65 
2563. 47 
2436. 20 
2304. 42 
2167. 70 

2025. 50 
1877. 18 
1721 . 95 
1558. 85 

1386. 61 
1203. 46 
1006 . 73 

791. 58 
546. 41 
211.63 
0 . 00 



-71.15 
-71.13 
-71.98 
-72.87 
-73.79 
-74.74 
-75.73 
-76.76 
-77.33 
-78.94 
-80.09 
-81 .30 
-82.56 
-83.87 
-85.24 
- 86.6 8 
-88.19 
-89 .78 
-91.44 
-93.19 
-95.04 
-97.00 
-99.06 
-101 .25 
-103.57 
-106.05 
-108.68 
-111.50 
-1 14.51 
-1 17.74 
- 121.21 
-124.94 
-128.98 
-133.33 
-138.05 
-143.15 
-148.67 
-154.64 
-161.05 
-167.88 
-175.05 
-182.41 
-1 89.67 
-196.50 
- 200.00 
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TABLE A-1 (CONT.) 



X X-VEL Y Y-VEL 



1011. 20 


10. 57 


1 C 21 .77 


25.14 


1046. 90 


35 . 8 2 


1082.72 


48 . 8 9 


1 131 . 6 1 


61.85 


1 153. 46 


74 . 4 9 


1267. 55 


86. 3 6 


1354.31 


97 .25 


1451 . 56 


107. 04 


1558 . 6 C 


115 . 75 


1674. 35 


123 . 4 3 


1797.78 


130.16 


1527.94 


136 . 06 


2 C 64.01 


141 . 24 


2205. 24 


145. 78 


2351 .02 


149 .78 


2500. 80 


153 . 3 1 


2654.1 1 


156.45 


2810. 56 


159 . 24 


2969 . 8 C 


161 .73 


3131. 52 


163 . 97 


3295 .45 


165 . 9 8 


3461 . 47 


167. 80 


3629. 27 


169 .46 


3798. 73 


170. 96 


3569.65 


172 .34 


4 14 2 . 0 3 


173. 60 


4315.63 


174. 76 


4490. 38 


175. 8 2 


4666. 21 


176 . 8 1 


4843. 02 


177 . 73 


5020.74 


178.57 


5199. 32 


179 . 36 


5378. 68 


180 . 1 0 


5558. 78 


180. 79 


5739 .57 


181 .44 


5 521. 0 1 


182. 04 


6103.05 


182 . 61 


6285.66 


183 . 14 


6468. 80 


183 . 65 


6652.45 


184 . 13 


6836.58 


184 . 58 


7021 . 16 


185. 0 0 


7206. 16 


185.41 


7391 . 57 


185. 79 



- 211. 63 


- 199.72 


- 295. 07 


- 198.41 


- 433 . 1 3 


- 196.77 


- 575. 23 


- 193.93 


-7254 57 


- 190.20 


- 879. 69 


- 1 85.61 


- 1035 . 28 


- 1 80.39 


- 1190. 48 


- 174.77 


- 1343 . 96 


- 168.94 


- 1494. 80 


- 163.10 


- 1642. 38 


- 157.37 


- 1786. 37 


- 1 51 . 35 


- 1926. 60 


- 146.58 


- 2063. 01 


- 141.61 


- 2195. 67 


- 136.92 


- 2324. 67 


- 132.54 


- 2450 . 14 


- 128.43 


- 2572. 25 


- 1 24.60 


- 2691 . 14 


- 121.01 


- 2806. 99 


- 117.66 


- 2919. 95 


- 114.52 


- 3030 . 17 


- 111.58 


- 3137. 81 


- 108.32 


- 3243. 01 


- 106.23 


- 3345 . 88 


- 103.79 


- 3446. 56 


- 101. 49 


- 3545 . 15 


- 99.32 


- 3641 . 77 


- 97.26 


- 3736. 51 


- 95.32 


- 3829. 47 


- 93.48 


- 3920. 72 


- 91.73 


- 4010. 36 


- 90.06 


- 4098. 45 


- 88.48 


- 41 85. 06 


- 86.97 


- 4270 . 26 


- 85.53 


- 4354 . 1 1 


- 84.15 


- 4436. 67 


- 82.83 


- 4517. 99 


- 81 .57 


- 4598 . 1 1 


- 80.36 


- 4677 . 09 - 


- 79.20 


- 4754. 98 


- 78.09 


- 4831 . 80 


- 77.01 


- 4907. 61 


- 75.98 


- 4982. 43 


- 74.99 


- 5056. 31 


- 74.03 
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Figure A-2 Track 4 



So each second: 



Tn — = *05 radians will be traversed 

lUTT 



Using the trigonometric identity: 
sin^(A) + cos^(A) = 1 
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And the equation for a circle: 



2 2 9 

X + y = c“ 

It follows that the arc of Figure A-2 is described by: 

x(k) = 4000 + lOOOsinC . 05k) 
y(k) = 4000 + lOOOcosC .05k) 

where the angle argument is in radians and k=0 , 1 , 2 , . . . 31 

seconds. The velocities v and v can be obtained from: 

X y 

V (k) = 50cos(.05k) 

X 

V (k) = -50sin(.05k) 

y 

using the same angle argument. 

The track values for the track arc are contained in Table 

A-2. 
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TABLE A-2 



NUMERICAL VALUES CF 



K 


X 


X- 


12 


3049 .9 S 


49 


13 


3099.96 


49 


14 


314 9.66 


49 


15 


3199.67 


49 


16 


3249.35 


49 


17 


3298. 88 


49 


18 


3348 . 22 


49 


IS 


3397.34 


49 


20 


3446.21 


48 


21 


3494.61 


48 


22 


354 3. C 9 


48 


23 


3591 . C4 


47 


24 


3638.62 


47 


25 


3685. eC 


46 


26 


3732.55 


46 


27 


3778.64 


46 


28 


3824 .64 


45 


29 


3869. S3 


45 


30 


3914.68 


44 


31 


3958.65 


43 


32 


4002.43 


43 


33 


4045.37 


42 


34 


4087 .67 


41 


35 


4129 .26 


41 


36 


4170. 19 


40 


37 


4210.37 


39 


38 


4249.79 


39 


“C 


4288 . 44 


38 


40 


4326 . 27 


37 


41 


4363.28 


36 


42 


4399 . 43 


35 


43 


4434 . 71 


34 


44 


4469. 10 


33 


45 


4502 .56 


33 


46 


4535. C9 


32 


47 


4566 . 65 


31 



SS EOR TRACK FOUR 



Y 


Y-V EL 


4999. 38 


-1.25 


4997. 50 


-2. 50 


4994. 38 


-3. 75 


4990. 01 


-4. 99 


4984. 40 


-6.23 


4977. 54 


-7.47 


4969. 45 


-8.71 


4960. 13 


-9. 93 


4949. 59 


-11.16 


4937. 82 


-1 2. 37 


4924. 65 


-1 3.53 


49 10. 67 


-14. 73 


4895. 30 


-1 5. 97 


4878. 75 


-17. 14 


4861, 02 


-18. 31 


4842. 12 


-19. 47 


4822. C8 


-20.62 


4800. 89 


-21. 75 


4778. 59 


-22. 87 


4755. 17 


-23. 97 


4730. 65 


-25.06 


4705 . 05 


-26. 13 


4678. 38 


-27. 19 


4650. 67 


-28. 23 


4621. S3 


-29.25 


4592. 17 


-30.26 


4561. 41 


-31.24 


4529. 68 


-32. 21 


4497. 00 


-33. 16 


4463. 38 


-3 4.0 8 


4428. 84 


-34.99 


4393. 41 


-35. 87 


4357. 1 1 


-36. 73 


4319. 97 


-37. 56 


4281. S9 


-38.38 


4243. 22 


-39. 17 



ST AT 

VEL 

. 98 

.94 

.86 

. 75 

. 6 1 

.44 

. 24 

.00 

.74 

.45 

. 1 2 

.77 

.38 

. 97 

. 53 

.05 

. 55 

.02 

. 46 

.88 

. 27 

. 63 

.96 

.27 

. 55 

.80 

.04 

.24 

. 42 

.58 

. 72 

.84 

. 93 

. 00 

. 05 

.08 
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U8 

49 

50 

51 

52 

53 

54 

c c 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

63 

69 

70 

71 

72 

73 

74 



TAEL5 A-2 (CONI,) 



X X-VEL Y Y-VEL 



4597. 24 


30,09 


4626. 63 


29 . 08 


4655. 4C 


28. 06 


4682 . 94 


27 . 02 


4709.43 


25. 9 5 


4734.65 


24.88 


4759. 18 


23. 79 


4782 . 4 1 


22 .68 


4804 . 54 


21 . 56 


4825.53 


20. 4 2 


4845. 38 


19. 2 8 


4864. C8 


18. 1 2 


4881.61 


16. 95 


4897 . 97 


15.77 


4913. 14 


14. 58 


4927 . 12 


13 .37 


4939.69 


12. 17 


4951 . 45 


10.95 


4961.79 


9.73 


4970. 9C 


8 .50 


4978.78 


7.26 


4985.43 


6 . 03 


4990.63 


4.78 


4994 .99 


3 .54 


4997.90 


2. 29 


4999.57 


1 . 04 


4999.98 


-.2 1 



4203. 67 


- 39.93 


4163. 37 


-40.67 


4122. 34 


-41 .39 


4080. 60 


-42.07 


4038. 20 


-42.74 


3995. 14 


-43.37 


3951. 46 


-43.98 


3907. 19 


-44.56 


3862. 35 


-45.11 


3816. 97 


-45.64 


3771. 09 


-46.13 


3724. 72 


-46.60 


3677. 69 


-47.04 


3630. 64 


-47.45 


3583. CO 


-47.83 


3535. 00 


-48.18 


3486. 66 


-48.50 


3438. 01 


-48.79 


3389. 10 


-49.04 


3339. 93 


-49.27 


3290. 56 


-49.47 


3241 . 01 


-49.64 


3191. 30 


-49.77 


3 141 . 47 


-49.87 


3091. 56 


-49.95 


3041 . 59 


-49.99 


2991. 59 


-50 .00 
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APPENDIX B 



COMPUTER PROGRAM EXPLANATION 



1. LINKAL 

The LINKAL program computes the filter gains, GAIN ( 4 , 2 ) , 
for the 4-state system and stores the gains in "LINGAIN 
.STORAG". The theoretical covariance of error matrix, PKK(4,4), 
is also computed and stored in "LINCOV . STORAG" . Several 
different sets of gain and covariance values were computed 
and stored for various values of measurement noise, RMAT ( 2 , 2 ) , 
and random forcing noise, COVW(2,2). 

2. LINEST 

The LINEST program retrieves the appropriate gain 
schedule from storage and computes the optimal estimate, 
XHAT(4,1), and the optimal one-step prediction, XHK1K(4,1). 

The following capabilities are contained in the LINEST program: 

a . Gating Scheme 

If the absolute difference (DIFF) between the one-step 
prediction, XHKIK and the noisy track value, ZMAT , is greater 
than the three-sigma gate, then the GAIN matrix is disregarded 
and XHAT is set equal to XHKIK. 

b . Track Noise 

By setting NOITRAK = 1, the random number generator, 

RND, and the resulting simulated noise produced, VI and V2, 
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are bypassed and the track values corrupted with noise, ZMAT , 
are retrieved from data file "NOITRAK . STORAG" . If hOITRA.K = l, 
then the random number generator is "reseeded" for each 
program run, producing a different set of noise values 
resulting in a unique noise corrupted track for each run. 

c . Statistics Window 

When WINDOW is set to 0, the filter state statistics 
are computed after each Monte Carlo run. The error mean, 
variance and position covariance are computed. If WINDOW 
is set to an integer, I, such that max k>I>0, the statistics 
will be computed based on the data compiled during the last 
I iterations of the simulation. If, for example, WINDOW=10, 
the computation of statistics will be based on the data 
obtained during the last ten iterations of k, and all 
previous data is disregarded. 

d . Error Nomalization 

By setting NORM = 1, the error (ERR), which is 
defined as the difference between the true track value (TRAK) 
and the estimate (XHAT), is normalized. For all other values 
of NORM the relative error is used in computing the 
statistics . 

e . Reinitializing Gains 

The program has the option of reinitializing the gains 
to G(0). This can be done by setting HIGH to an integer I, 
such that max k>I>0. If the surface area of the ellipse 
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increases I consecutive times, indicating filter divergence, 
the gains are reinitialized to G(0), If reinitializing 
gains is not desired as an option, then HIGH is set to some 
arbitrary large number greater than max k. 

3 . EXTKF 

The EXTKF program computes the optimal estimate, XHAT , 
and the optimal one-step prediction, XHKIK, for the 4-state 
nonlinear tracking problem. The EXTKF program has the 
options: gating scheme, track noise, statistics window, 

and error normalization as described for the FINEST program. 
EXTKF has the following additional options : 

a . Reinitializing the Filter 

The program has the option of reinitializing the 

covariance of one-step prediction error matrix, PKIK, by 

setting it to lO^xP^, where n is zero or some small integer. 

If HIGH is set to an integer I, such that max k>I>0, and the 

surface area of the ellipse increases I consecutive times, 

then the PKIK matrix will be reset to lO^xP . 

— o 

b . Adaptive Q 

EXTKF has the option of automatically increasing or 
decreasing the state excitation matrix, Q_, under certain 
conditions by changing the value of the covariance of 
excitation noise, COVWo INCREASE and DECREASE are set to 
integers I and J respectively, such that max k>I,J>0. If 
the surface area of the error ellipse increases I 
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consecutive times, COVW is doubled or increased by some 
factor. Conversely, if the surface area of the error ellipse 
decreases J consecutive times, COVW is halved or decreased 
by some factor. If COVW is changed, increased for example, 
the value of COVW can be changed again later in the run, 
if the criteria for either increasing or decreasing is met. 

If the adaptive Q is not desired then INCREASE and DECREASE 
are set to some large number greater than k. 
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